% fig3 李指数图(LEs)

clear; clc; close all;

%% —— 参数设置 —— 
g11 = 1; g12 = 1; g21 = -2; g22 = 2;   % 连接权重
mu_vec = linspace(0,1,500);           % μ 从 0 到 1，共 500 个采样点
nMu    = numel(mu_vec);

% 初始条件
x0 = 0.1;  y0 = 0.1;

% 迭代总步数与热身步数
Ntot   = 3000;   % 总迭代步
Ntrans = 1000;   % 去除热身

% 预分配
LE1 = zeros(nMu,1);
LE2 = zeros(nMu,1);

%% —— 并行计算李雅普诺夫指数 —— 
parfor idx = 1:nMu
    mu = mu_vec(idx);

    % 初始化相空间与两条正交扰动向量
    x = x0; y = y0;
    v1 = [1;0];   v2 = [0;1];
    s1 = 0; s2 = 0;

    for k = 1:Ntot
        % —— 离散映射迭代 —— 
        x_next = mu*x + g11*sin(x) + g12*sin(y);
        y_next = mu*y + g21*sin(x) + g22*sin(y);

        % —— 雅可比矩阵 —— 
        J = [mu + g11*cos(x),      g12*cos(y);
             g21*cos(x),      mu + g22*cos(y)];

        % —— 同步扰动迭代 —— 
        w1 = J * v1;
        w2 = J * v2;

        % —— 正交分解（Gram–Schmidt） —— 
        % 第一向量归一
        norm1 = norm(w1);
        v1 = w1 / norm1;
        % 第二向量去除 v1 分量再归一
        w2 = w2 - (v1' * w2) * v1;
        norm2 = norm(w2);
        v2 = w2 / norm2;

        % —— 累加对数伸缩 —— 
        if k > Ntrans
            s1 = s1 + log(norm1);
            s2 = s2 + log(norm2);
        end

        % 准备下一步
        x = x_next;  y = y_next;
    end

    % 平均并存储 LE
    M = Ntot - Ntrans;
    LE1(idx) = s1 / M;
    LE2(idx) = s2 / M;
end

%% —— 绘制 LE 光谱 —— 
figure;
plot(mu_vec, LE1, 'r-', 'LineWidth', 1.2); hold on;
plot(mu_vec, LE2, 'b-', 'LineWidth', 1.2);
xlabel('\mu', 'FontSize',14);
ylabel('LEs',  'FontSize',14);
title('第3图(b)：\mu–LE 光谱', 'FontSize',16);
legend('LE_1','LE_2','Location','northeast');
xlim([0,1]);
grid on;
